###################################################################################
# R FUNCTIONS 
# "When do the Wealthy Support Redistribution?," BJPS, forthcoming 
# G. Feierherd, L. Schiumerini, and S. Stokes 
###################################################################################

# model_iv is necessary to replicate figures 3 and 4 in the main text. 
# This function runs the CACE or instrumental variable model with interactions. 
# It shows the mean for each group and the marginal effects with its standard errors.
# It allows to specify the data, dependent variable, experimental treatment, the baseline, policy treatment, and the IV.

model_iv <- function(data=data, dep, treat, base, lost, iv) {
  
  # creating temporal DS
  temp <- data[treat ==1 | base ==1, ]
  dep.temp <- dep[treat ==1 | base ==1]
  treat.temp  <- treat[treat ==1 | base ==1]
  lost.temp  <- lost[treat ==1 | base ==1]
  iv.temp  <- iv[treat ==1 | base ==1]
  
  # Vs. Base
  cace <- ivreg(dep.temp~treat.temp+lost.temp+treat.temp*lost.temp,
                ~treat.temp+iv.temp+treat.temp*iv.temp, data=temp)
  cace.SE <- vcovHC(cace, type = "HC3") 
  b_11 <- cace$coefficients[1] 
  b_12 <- b_11 + cace$coefficients[2] 
  b_21 <- b_11 + cace$coefficients[3] 
  b_22 <- b_11 + cace$coefficients[2] + cace$coefficients[3] + cace$coefficients[4] 
  Marg1 <- cace$coefficients[2] 
  Marg2 <- cace$coefficients[2] + cace$coefficients[4] 
  se_11  <- sqrt(cace.SE[1,1])
  se_12  <- sqrt(cace.SE[1,1] + cace.SE[2,2] + 2*cace.SE[1,2])
  se_21  <- sqrt(cace.SE[1,1] + cace.SE[3,3] + 2*cace.SE[1,3])
  se_22  <- sqrt(cace.SE[1,1] + cace.SE[2,2] + cace.SE[3,3] + cace.SE[4,4] + 2*cace.SE[1,2] + 2*cace.SE[1,3] +  2*cace.SE[1,4]  + 2*cace.SE[2,3] + 2*cace.SE[2,4] + 2*cace.SE[3,4])
  se_Marg1 <- sqrt(cace.SE[3,3])
  se_Marg2 <- sqrt(cace.SE[3,3] + cace.SE[4, 4] + 2*cace.SE[3, 4])
  #dm.n <- summary(dm)$df[2]
  dm1 <- c(b_11,b_12,Marg1)
  dm2 <- c(b_21,b_22,Marg2)
  se1 <- c(se_11,se_12,se_Marg1)
  se2 <- c(se_21,se_22,se_Marg2)
  results <- as.data.frame(round(rbind(dm1, se1, dm2, se2), digits=3))
  colnames(results) <- c("Control", "Treat", "Marginal Effect")
  rownames(results) <- c("Kept", "SE1","Lost", "SE2")
  
  return(results)
  
}

# group_iv is necessary to replicate figures A6 and A7 in the online appendix. 
# It plots the outcome of model_iv

groups.iv <- function(data,treat,base,lost,iv,indep.lab="",main.title="",text=1,a=-.75,b=.75) {
  
  # Upper and lower confidence bounds
  
  z_score = qnorm(1 - ((1 - .95)/2))
  
  pol <- model_iv(data,data$p1_1z,treat,data$placebo,data$lost,data$nat)
  
  pu21 <- pol[3,1] + z_score*pol[4,1]
  pu22 <- pol[3,2] + z_score*pol[4,2]
  pl21 <- pol[3,1] - z_score*pol[4,1]
  pl22 <- pol[3,2] - z_score*pol[4,2]
  
  red <- model_iv(data,data$p2_1z,treat,data$placebo,data$lost,data$nat)
  
  ru21 <- red[3,1] + z_score*red[4,1]
  ru22 <- red[3,2] + z_score*red[4,2]
  rl21 <- red[3,1] - z_score*red[4,1]
  rl22 <- red[3,2] - z_score*red[4,2]
  
  une <- model_iv(data,data$p2_2z,treat,data$placebo,data$lost,data$nat)
  
  uu21 <- une[3,1] + z_score*une[4,1]
  uu22 <- une[3,2] + z_score*une[4,2]
  ul21 <- une[3,1] - z_score*une[4,1]
  ul22 <- une[3,2] - z_score*une[4,2]
  
  # Initialize plotting window  
  plot(x=c(), y=c(), ylim=c(-1.25, 1.25), xlim=c(0,12), xlab="", 
       ylab="Standardized mean",  xaxt="n", type = "n", bty='n', axes=T, main=main.title)
  
  abline(h=0, lty=3)
  
  axis(1,at=c(2,6,10),labels=c("Opinion about\n the policy","Redistribution",
                               "Unemployment\ninsurance"), mgp=c(3, 2, 0))
  # Price Hike 
  
  points(x=1, y=pol[3,1], pch=16, col="red", cex=1.5)
  points(x=3, y=pol[3,2], pch=15, col="blue", cex=1.5)
  
  lines(x=c(1,1), y=c(pl21, pu21), col="red", lty=1, lwd=3)
  lines(x=c(3,3), y=c(pl22, pu22), lty=1, col="blue", lwd=3)
  
  segments(1.05,pol[3,1],2.95,pol[3,2], lwd=4, lty=3, col="gray")
  
  text(1-.5, pol[3,1], round(pol[3,1],2),srt=0, cex=1)
  text(3+.5, pol[3,2], round(pol[3,2],2),srt=0, cex=1)
  
  sig(ate=pol[3,3],se=pol[4,3],x=1,y=1,f=1,est="CATE = ")
  
  # Redistribution
  
  points(x=5, y=red[3,1], pch=16, col="red", cex=1.5)
  points(x=7, y=red[3,2], pch=15, col="blue", cex=1.5)
  
  lines(x=c(5,5), y=c(rl21, ru21), col="red", lty=1, lwd=3)
  lines(x=c(7,7), y=c(rl22, ru22), lty=1, col="blue", lwd=3)
  
  segments(5.05,red[3,1],6.95,red[3,2], lwd=4, lty=3, col="gray")
  
  text(5-.5, red[3,1],round(red[3,1],2),srt=0, cex=1)
  text(7+.5, red[3,2],round(red[3,2],2),srt=0, cex=1)
  
  sig(ate=red[3,3],se=red[4,3],x=6,y=1,f=1,est="CATE = ")
  
  # Unemployment insurance
  
  points(x=9, y=une[3,1], pch=16, col="red", cex=1.5)
  points(x=11, y=une[3,2], pch=15, col="blue", cex=1.5)
  
  lines(x=c(9,9), y=c(ul21, uu21), col="red", lty=1, lwd=3)
  lines(x=c(11,11), y=c(ul22, uu22), lty=1, col="blue", lwd=3)
  
  segments(9.05,une[3,1],10.95,une[3,2], lwd=4, lty=3, col="gray")
  
  text(9-.5, une[3,1],round(une[3,1],2),srt=0, cex=1)
  text(11+.5, une[3,2],round(une[3,2],2),srt=0, cex=1)
  
  sig(ate=une[3,3],se=une[4,3],x=10,y=1,f=1,est="CATE = ")
  
  
}

# sig adds stars indicating statitical significance in plots and tables 

sig<-function(ate,se,x,y,f=.5,est){
  
  if (2*pnorm(-abs(ate/se))<.001){
    
    
    text(x, y, paste(est,round(ate,digits=2)," (",round(se,digits=2),")","***",sep=""),cex=f,font=4)
    
  }else{
    
    
    if (2*pnorm(-abs(ate/se))<.05){
      
      
      text(x, y , paste(est,round(ate,digits=2)," (",round(se,digits=2),")","**",sep=""),cex=f,font=4)
      
    }else{
      
      if (2*pnorm(-abs(ate/se))<.1){
        
        
        text(x, y , paste(est,round(ate,digits=2)," (",round(se,digits=2),")","*",sep=""),cex=f,font=4)
        
      }
      
      else{
        
        
        
        
        text(x, y , paste(est,round(ate,digits=2)," (",round(se,digits=2),")",sep=""),cex=f,font=4)
        
        
      }
    }
  }
  
  
}

# pmean calculates row means 

pmean <- function(..., na.rm=FALSE) { 
  x <- list(...)
  rowMeans(matrix(unlist(x), ncol=length(x)), na.rm=na.rm)
}

# bal.surv calculates balance statistics for survey experimental treatments 

bal.surv <- function(outcome.var, base.var, treat1.var, treat2.var, data){
  
  library(car)
  data$outcomevar <- data[,outcome.var]
  data$base <- data[,base.var]
  
  data$treat1 <- data[,treat1.var]
  data1 <- data[data$base==1 | data$treat1==1,]
  
  data$treat2 <- data[,treat2.var]
  data2 <- data[data$base==1 | data$treat2==1,]
  
  dm.results1 <- lm(outcomevar ~treat1,data = data1)
  dm.out1 <-  c(coef(dm.results1)[[2]], sqrt(vcovHC(dm.results1, type = "HC2")[2,2]), dm.results1$df.residual + 2, sd(data1$outcomevar, na.rm = TRUE), coef(dm.results1)[[2]]+coef(dm.results1)[[1]], coef(dm.results1)[[1]])
  names(dm.out1) <- c("Estimate", "Standard Error", "N", "SD", "T", "C") 
  
  dm.results2 <- lm(outcomevar ~treat2,data = data2)
  dm.out2 <-  c(coef(dm.results2)[[2]], sqrt(vcovHC(dm.results2, type = "HC2")[2,2]), dm.results2$df.residual + 2, sd(data2$outcomevar, na.rm = TRUE), coef(dm.results2)[[2]]+coef(dm.results2)[[1]], coef(dm.results2)[[1]])
  names(dm.out2) <- c("Estimate", "Standard Error", "N", "SD", "T", "C") 
  
  list(dm.results1 = dm.out1,dm.results2 = dm.out2)
}

# model_itt is necessary to replicate Figure A5: ITT analysis of Resentment (H2) and Empathy (H3)

model_itt <- function(data=data, dep, treat, base, lost) {
  
  # creating temporary DS
  temp <- data[treat ==1 | base ==1, ]
  dep.temp <- dep[treat ==1 | base ==1]
  treat.temp  <- treat[treat ==1 | base ==1]
  lost.temp  <- lost[treat ==1 | base ==1]
  
  # Vs. Base
  cace <- lm(dep.temp~treat.temp+lost.temp+treat.temp*lost.temp, data=temp)
  cace.SE <- vcovHC(cace, type = "HC2") 
  b_11 <- cace$coefficients[1] 
  b_12 <- b_11 + cace$coefficients[2] 
  b_21 <- b_11 + cace$coefficients[3] 
  b_22 <- b_11 + cace$coefficients[2] + cace$coefficients[3] + cace$coefficients[4] 
  Marg1 <- cace$coefficients[2] 
  Marg2 <- cace$coefficients[2] + cace$coefficients[4] 
  se_11  <- sqrt(cace.SE[1,1])
  se_12  <- sqrt(cace.SE[1,1] + cace.SE[2,2] + 2*cace.SE[1,2])
  se_21  <- sqrt(cace.SE[1,1] + cace.SE[3,3] + 2*cace.SE[1,3])
  se_22  <- sqrt(cace.SE[1,1] + cace.SE[2,2] + cace.SE[3,3] + cace.SE[4,4] + 2*cace.SE[1,2] + 2*cace.SE[1,3] +  2*cace.SE[1,4]  + 2*cace.SE[2,3] + 2*cace.SE[2,4] + 2*cace.SE[3,4])
  se_Marg1 <- sqrt(cace.SE[3,3])
  se_Marg2 <- sqrt(cace.SE[3,3] + cace.SE[4, 4] + 2*cace.SE[3, 4])
  dm1 <- c(b_11,b_12,Marg1)
  dm2 <- c(b_21,b_22,Marg2)
  se1 <- c(se_11,se_12,se_Marg1)
  se2 <- c(se_21,se_22,se_Marg2)
  results <- as.data.frame(round(rbind(dm1, se1, dm2, se2), digits=3))
  colnames(results) <- c("Control", "Treat", "Marginal Effect")
  rownames(results) <- c("Kept", "SE1","Lost", "SE2")
  
  return(results)
  
}

